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This  paper  studies  the  applicability  of  the 
path  integral  solution  technique  for  esti- 
mating extreme  response  of  nonlinear  dy- 
namic oscillators  whose  equations  of  mo- 
tion can  be  modelled  by  the  use  of  Ito 
stochastic  differential  equations.  The  state 
vector  process  associated  with  such  a 
model  is  generally  a diffusion  process, 
and  the  probability  density  function  of  the 
state  vector  thus  satisfies  the  Fokker- 
Pl an ck- Kolmogorov  equation.  It  is  shown 
that  the  path  integral  solution  technique 
combined  with  an  appropriate  numerical 
scheme  constitutes  a powerful  method  for 
solving  the  Fokker-Planck-Kolmogorov 
equation  with  natural  boundary  condi- 


tions. With  the  calculated  probability  den- 
sity function  of  the  state  vector  in  hand, 
one  can  proceed  to  calculate  the  required 
quantities  for  estimating  extreme  re- 
sponse. The  proposed  method  distinguishes 
itself  by  remarkably  high  accuracy  and 
numerical  robustness.  These  features  are 
highlighted  by  application  to  example 
studies  of  nonlinear  oscillators  excited  by 
white  noise. 
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1.  Introduction 

An  important  element  in  the  safety  assessment  of 
many  engineering  systems,  is  the  task  of  estimating  the 
probability  of  extreme  events  that  may  jeopardize  the 
structure  in  some  specified  sense.  Very  often,  this  prob- 
lem can  be  formulated  as  finding  the  probability  that 
some  time  varying  random  quantity  does  not  exceed  a 
specified  capacity  level  during  a given  time  period. 
Stated  this  way,  the  problem  typically  reduces  to  a study 
of  the  extreme  values  of  a stochastic  process  originating 
as  the  response  of  a system  subjected  to  some  stochastic 
loading  process. 

In  this  paper  the  focus  will  be  on  the  problem  of 
estimating  the  extreme  response  of  nonlinear  dynamic 
systems  subjected  to  random  forcing  processes.  In  re- 
cent years  the  methods  of  time  domain  Monte  Carlo 
simulations,  see,  e.g..  Refs.  [1-5],  have  received  consid- 


erable attention  as  a tool  for  estimating  response  statis- 
tics. These  methods  are  versatile  and  attractive  in  the 
sense  that  nonlinearities  can  be  easily  dealt  with.  The 
main  drawback  at  present  is  the  large  CPU  times  needed 
for  accurate  prediction  of  extreme  responses.  Even  if 
this  issue  seems  to  become  less  of  an  obstacle  every 
year,  portending  perhaps  that  such  methods  may  domi- 
nate practical  estimation  of  response  statistics  of  nonlin- 
ear systems  in  the  not  too  distant  future,  it  will  still  be 
desirable  to  have  available  alternative  methods  of  calcu- 
lating the  response  statistics,  both  simplified  and  more 
elaborate.  Here  we  shall  explore  a method  based  on  the 
theory  of  Markov  diffusion  processes.  The  justification 
for  using  this  theory  is  related  to  the  fact  that  the  re- 
sponse of  nonlinear  dynamic  systems  to  broad  band 
random  excitation  can  very  often  be  accurately  de- 
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scribed  by  applying  the  theory  of  multidimensional 
Markov  processes.  By  this,  the  extensive  theory  of 
Markov  diffusion  processes  can  be  brought  to  bear  on 
these  problems.  In  particular,  it  can  be  shown  that  the 
probability  law  of  response  quantities  can  be  derived  by 
solving  a partial  differential  equation,  viz.,  the  Fokker- 
Planck  (-Kolmogorov)  (FPK)  equation,  see  Refs.  [6,7]. 
In  most  cases  of  practical  interest,  this  equation  has  to  be 
solved  numerically. 

In  the  next  section  we  shall  describe  a method  for 
solving  the  FPK  equation  that  is  based  on  a formal 
solution  of  the  same  equation.  This  solution  is  obtained 
by  invoking  the  fact  that  a Markov  diffusion  process 
locally  looks  like  a Brownian  motion.  By  using  the 
Markov  property,  the  global  solution  can  then  be  con- 
structed by  linking  the  local  solutions,  which  are  known 
explicitly.  The  obtained  solution  is  generally  known  as  a 
path  integral  solution  (PIS).  The  reader  is  referred  to 
Ref.  [7]  for  a further  discussion.  One  of  the  first  efforts 
to  exploit  the  PIS  method  explicitly  in  developing  nu- 
merical solution  algorithms  is  described  in  Ref.  [8]. 
Subsequently,  other  authors  have  also  used  the  PIS  ap- 
proach to  solve  various  random  vibration  problems,  cf. 
Refs.  [9-13]. 

Before  embarking  on  a description  of  the  PIS  method, 
it  is  expedient  to  briefly  show  how  the  obtained  solutions 
are  used  in  an  extreme  value  analysis.  Assuming  that  the 
response  quantity  of  interest  is  a scalar  (real)  stationary 
stochastic  process,  Z{t)  say,  the  PIS  method  typically 
provides  a numerical  estimate  of  the  joint  probability 
density  function  (PDF) /zz(*,*)  of  Z{t)  and  Z{t)=dZ(t)l 
dt . It  is  now  assumed  that  the  mean  level  upcrossing  rate 
vzi*)  of  Z(r)  can  be  calculated  from  Rice’s  formula  as 
follows 

vHz)=jyfzz{z,y)dy.  (1) 

0 

Adopting  the  assumption  that  upcrossing  of  high  levels 
are  statistically  independent  events,  which  leads  to  Pois- 
son distributed  crossings,  it  follows  that  an  asymptotic 
approximation  of  the  probability  distribution  function  of 
the  extreme  value  of  the  process  Z{t)  during  a time  T, 
denoted  by  MiT)  (=sup{Z(0;0:£t:£7’}),  is  given  by 

Prob{M{T)<z}=e\p{-Vz{z)T}  (T-^oc).  (2) 

The  accuracy  of  Eq.  (2)  depends  to  a large  extent  on 
the  effective  bandwidth  of  the  response  process  Z(/). 
Decreasing  bandwidth  leads  eventually  to  a significant 
clumping  effect  of  large  response  peaks,  invalidating 
the  assumption  of  statistically  independent  upcrossing  of 
high  levels.  Methods  that  aim  at  correcting  for  this  effect 
have  been  proposed  for  Gaussian  (Refs.  [14,15])  and 


non-Gaussian  (Ref.  [16])  processes.  However,  this  point 
will  not  be  pursued  any  further  here.  We  shall  assume 
that  Eq.  (2)  provides  an  acceptable  approximation. 
Hence,  the  central  parameter  to  be  determined  is  the 
upcrossing  frequency  i4(*),  which  is  easily  calculated 
once  the  joint  PDF /zz(*/)  of  Z(t)  and  Z{t)=dZ{t)ldt  has 
been  made  available.  In  the  next  section  it  is  shown  how 
this  PDF  can  be  calculated  for  the  response  of  a wide 
range  of  nonlinear  oscillators  subjected  to  white  noise 
or  filtered  white  noise  loading. 

2.  The  Path  Integral  Solution 

The  path  integral  solution  (PIS)  method  is  suitable  for  > 
calculating  the  joint  probability  density  function  (PDF) 
of  a vector  process  2f(r)=[Xi(r),...,X„(r)]^  (T-transposi- 
tion)  satisfying  a stochastic  differential  equation  of  the 
following  form,  cf.  Ref.  [6], 

dXit)=m[Xit)]dt+Q[X{t)]dWit).  (3)  ! 

Here  /w  (•)=[/« i(»),...,/n„(»)]^,  nij{*)  denotes  a real  func- 
tion of  n real  variables.  Qi*)=(qij{*))  denotes  an  nXm- 
matrix  where  each  is  a real  function  of  n real 

variables.  W(r)=[VFi(r),...,VT„(f)]^  where  Wjit), 
y=l,...,m  are  standard,  real  Brownian  motion  processes, 
which  are  mutually  independent,  see  e.g..  Refs.  [6,7]. 

That  is,  £'[VVj(f)]=0  and 

£:[dW,(r)dVT;(r-f-T)]=8y8,,,^.^dT,  /J=l,...,m,  (4) 

where  8^=1  for  x=y,  8tj,=0  for  X7^y.  Equation  (4) 
is  a short-hand  notation  for  the  relation 

/2(5,0dW,(5)dW;(r)]=8,J /2(r,/)dr,  where  /?(*,*) 

i 

is  a non-random  function.  i 

Equation  (3)  is  interpreted  here  as  an  Ito  stochastic 
differential  equation  (SDE).  Since  it  is  often  relevant  to 
consider  Eq.  (3)  as  being  obtained  as  a limit  of  equations 
with  band  limited  noise  processes,  it  may  happen  that 
m(«)  should  contain  correction  terms  to  ensure  a consis- 
tent limiting  solution,  cf.  Ref.  [6].  It  is  assumed  here  that 
this  consideration  has  already  been  made,  and  that  Eq. 

(3)  has  the  final  form  to  be  used  subsequently. 

It  is  demonstrated  in  Ref.  [6]  that  the  solution  A" (r)  to 
Eq.  (3)  is  a Markov  vector  process.  Its  transition  proba- 
bility density  function  (TPD),  p(x,t  is  defined 

by  the  equation 

Prob{X{t)  E.  A I A'(r')=jc'}=j'«»»j'p(jc,/ |x',t')dx,  (5)  | 

where  A C R"  is  some  event,  x,  x'  E R”,  djc=dxi...dr„.  i 
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Provided  that  m(»)  and  Q(»)  satisfy  certain  regularity 
conditions,  see  Ref.  [6],  it  can  be  proved  that  the  TPD 
p{x,t  I x' ,t')  {t^t'^0)  is  the  solution  of  a partial  differ- 
ential equation  of  the  form 


O denote  appropriate  zero-matrices  and  G(»)  denotes  an 
(«— r)X(«— r)-matrix  function  with  elements  gy(*), 
ij=r+l,...,n.  G(*)  will  be  called  the  reduced  diffusion 
matrix.  Equation  (6)  can  now  be  rewritten  as 


— p(x,r  lx',t')=-S  ■^[m,{x)p{x,t\x',t')] 


1 " " 
^ /-I  j-i 


(6) 


I x',t')=-S  -^[mi{x)p(x,t  I x',t')] 

+ i [gij{x)p{x,t\x',t')].  (12) 

^ i-r+\  j-r+l 


where  G(x)={gijix))=Q{x}Q{xf=i^'^_^qj,qji),  and 

with  initial  condition  pix,t'  | x',r')=8(x— x').  G(*)  will 
be  called  the  diffusion  matrix  and  Eq.  (6)  will  be  re- 
ferred to  as  the  Fokker-Planck-Kolmogorov  (FPK) 
equation.  Since  clearly  Prot>  {X(t)  E R''fx(f')=x' }=l 
the  TPD  satisfies  the  following  normalization  condition 

•yjp(x,r  |x',r')dx=l.  (7) 

Let  f(x,t)  denote  the  PDF  of  the  random  vector  A'(t). 
If  /(x,f')=w(x)  for  some  initial  PDF  w(x),  then  it  is 
recognized  from  Eq.  (6)  and  the  relation 


Proceeding  in  a manner  similar  to  the  derivations  in 
Ref.  [7],  it  can  be  shown  that  the  TPD  for  small  values 
of  r{=t—t')  is  given  by  the  following  expression,  which 
is  correct  up  to  terms  of  order 

p{x,t+T  I x',/)={n5(x,-X/'-m,(x')T)} 

/-I 

. n n 

1 G(x')  I '2  ‘expl-—  X 2 

i^r+]  y-r-bl 

{x-x'-m, (x' )t)  [G(x' )■ ' l-rj-r{Xj-Xj-mj{x' )t) }, 

(13) 


x',r')w(x')dx' 


(8) 


that  f{x,t)  itself  is  a solution  of  Eq.  (6)  satisfying  the 
initial  condition  f{x,t')=w{x). 

In  this  paper  we  shall  be  interested  primarily  in  sta- 
tionary solutions  /j(x)  to  Eq.  (6),  that  is 

/,(x)=lim/(x,r)=lim  p(x,t\  x',t’}  (9) 

r^co  /_>oo 

provided  they  exist.  Even  when  both  limits  exist,  it  is 
clear  that  lim f{x,t)  provides  the  faster  convergence 
when  the  initial  condition /(x,/')~/i(x).  This  comment 
is  relevant  to  the  numerical  implementation  of  the  PIS 
method,  and  will  be  discussed  below. 

To  obtain  the  PIS  appropriate  for  the  dynamic  systems 
studied  in  this  paper,  it  is  necessary  to  be  more  specific 
on  the  structure  of  the  matrix  function  Q{*).  In  particu- 
lar, it  will  be  assumed  that  the  first  r rows  of  Q(*)  are 
zero,  that  is 


where  | G I denotes  the  determinant  of  the  reduced  dif- 
fusion matrix  G,  assumed  to  be  positive  definite.  This 
implies  that  | G | > 0.  [G“'],j  denotes  the  element  in  po- 
sition ij  of  the  inverse  matrix  of  G.  As  shown  in  Ref. 
[7],  the  expression  given  by  Eq.  (13)  is  not  unique,  but 
seems  to  be  well  suited  for  our  purpose. 

Having  obtained  an  explicit  expression  for  the  TPD 
for  a short  time  step,  one  can  now  invoke  the  Markov 
property.  This  allows  a TPD  over  a time  interval  of 
arbitrary  length  to  be  expressed  in  terms  of  a product  of 
short-time  TPDs.  By  dividing  a given  time  interval 
into  /V  small  time  intervals  of  length  T={t—t')/N,  it  is 
found  that  itj=t'+jT,  x=x'^,  t'=k),  x'=x''”) 


p(x,r  I x',r')=  f***  ff7p(x‘^\/^  Ix^  '’,/^-i)dx''’»»»dx‘'^ 

(14) 


Similarly,  with  an  initial  PDF/(x',r')=w(x'),  the  PDF 
f{x,t)  will  be  given  by 


q,j{»)^  for  {r<n)  (10) 

and  that  qij{*)^0  for  at  least  one  j for  every  i-r+],...,n. 
This  implies  that  the  diffusion  matrix  G(»)  assumes  the 
form 


G(*)= 


O 

O 


o 

G(*) 


(11) 


x‘^-'>,0_,)w(x<"')dx''’>— dx‘^-'>. 

(15) 


Hence,  by  combining  Eq.  ( 1 3)  with  Eqs.  ( 1 4)  or  ( 1 5), 
a formal  (approximate)  solution  of  the  FPK  equation  can 
be  written.  Equations  (14)  and  (15),  which  are  often 
referred  to  as  PIS,  constitute  the  core  of  the  numerical 
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solution  procedure  to  be  described  subsequently.  It  is 
realized  that  a numerical  solution  according  to  this 
method,  automatically  provides  the  evolution  in  time  of 
the  (conditional)  PDF  of  the  Markov  process  X(t)  from 
given  start  conditions  in  terms  of  an  initial  density 
fix'  ,t')=w{x'),  including  the  degenerate  case 
fix'  ,t')=hix'—Xo),  for  some  starting  point  Xq.  It  is  also 
worth  noting  how  the  PIS  relates  to  the  physics  of  the 
dynamic  model,  which  is  expressed  through  the  coeffi- 
cients /«;(•)  and  qtji*),  cf.  Eq.  (3).  The  evolution  in  time 
of  the  PDF  as  expressed  by  the  PIS,  is  seen  to  be  directly 
determined  by  these  coefficients  in  an  explicit  manner. 
This  fact  is  a very  important  advantage  of  the  PIS 
method,  and  reveals  its  fundamental  physical  signifi- 
cance. 

3.  Numerical  Implementation 

In  the  numerical  implementation,  the  PIS  is  obtained 
by  an  iteration  process  based  on  the  Chapman-Kol- 
mogorov  equation  expressed  as 

pix^^,tj  \x^~'\tj-i)pix^~'\ 

tj-Ax',t')(ix^-'\  (16) 

The  discretization  of  state  space  for  the  numerical 
solution  makes  it  appropriate  to  employ  an  interpolation 
and  smoothing  procedure  to  increase  the  numerical  effi- 
ciency. It  was  found  that  application  of  cubic  B-splines, 
as  detailed  in  Ref.  [17],  offered  the  desired  accuracy  and 
smoothness  for  the  type  of  problems  considered  in  this 
paper.  This  procedure  was  used  as  follows.  At  each  time 
step  pix^~'\tj-i  I x' ,t')  is  represented  as  a cubic 

B-spline  series  in  the  following  manner 

M|  M„ 

p ix^-'Ktj ! X' ,/' )=  2 •••X ) 

tl-l 

<^B,ix^-\  (17) 

(-1  ' 

where  Mj=number  of  grid  points  for  the  i’th  state  vari- 
able X,,  is  a tensor  product  basis  of  cubic 

B-splines  and  is  the  set  of  interpola- 

tion coefficients  associated  with  time  t;-,.  It  is  assumed 
that  each  set  {Bt,(*)}I^li,  /=1, •••,«,  is  a basis  of  cubic 
B-splines  associated  with  the  knot  sequence  determined 
by  the  grid  points  for  the  i’th  variable  x,.  The  tensor 
product  B-spline  is  defined  by 

n " 

<S)B,fx)=Y[B,,ix,).  (18) 

,-i 


The  representation  of  p(x^“’Vj-i  l^r'^F)  by  B-splines 
makes  it  possible  to  retain  high  numerical  accuracy  even 
with  a fairly  coarse  basic  grid  if  p(x^“‘’,(,-i  | x',t')  is  not 
too  singular.  By  substituting  from  Eq.  (17)  into  Eq.  (16), 
Eq.  (19)  is  obtained 

M|  M„ 

p(x^Vj  U',r')=2***E  r^~'\ki,»»\k„) 

*1-1  *n-l 

j J P 1 ' ’f- 1 ) ^ Bn,  ix^- ' ')dx^" ' (19) 

It  is  seen  from  Eq.  (13)  that  since  mji»)  and  g/ji*)  are 
not  functions  of  time  t,  the  TPDs  cannot  depend  on 
absolute  time,  but  only  on  the  time  increment.  Markov 
processes  whose  TPDs  have  this  property,  are  called 
homogeneous.  It  follows  that 

pix^\tj  I x^-'\tj-,)=p(x^\r  I x^-'>,0),;=l,2,...,  (20) 

which  holds  for  any  t—tj-\=f^Q. 

From  Eqs.  (19)  and  (20)  it  is  seen  that  for  a fixed 
value  of  the  time  increment  t,  each  of  the  integrals  on 
the  right  hand  side  of  Eq.  (19)  need  to  be  calculated  only 
once,  and  can  be  stored  for  repeated  use.  That  is,  the 
following  parameters  are  calculated  initially  and  stored 

f*  * * I ® 5,,(x^-'>)dx^-'>. 

J R"  J 1-1 

(21) 

Here,  the  index  /,  , /=1, ...,  n,  refers  to  grid  point  number 
/,  for  the  state  space  variable  x,  . It  may  be  noted  here  that 
due  to  the  properties  of  the  TPD  for  small  time  incre- 
ments T,  the  tensor  Bl\"f„  has  a strongly  banded  charac- 
ter with  the  elements  decreasing  rapidly  away  from  the 
main  diagonal  k\=l\,...,k„=h-  This  has  important  impli- 
cations for  the  efficiency  of  the  computer  program.  Let 
P^lj=P(xu].j„)f  ! x',t').  Then  Eq.  (19)  can  be  rewritten 
as 

M]  M„ 

I^-'\ku...X)Bl\X.  (22) 

it  1-1  k„~\ 

Having  calculated  the  TPD  p(x^\(,  | x',t')  at  the  grid 
points  by  using  Eq.  (22),  a spline  interpolation  is  again 
carried  out  and  a new  set  of  interpolation  coefficients 
{r^\k ,kn)}kri  are  calculated.  This  provides  an  up- 

dated representation  of  the  TPD  for  time  step  j,  cf.  Eq. 
(17).  For  each  time  step,  the  normalization  condition  Eq. 
(7)  is  checked.  That  is,  if 

r 

• • • p(x^V;  I x',/')dx^^=X  "‘S  J^n) 

J i,-i  VI 

n [ Bk,ix)6x-=qj  (23) 


468 


Volume  99,  Number  4,  July-August  1994 

Journal  of  Research  of  the  National  Institute  of  Standards  and  Technology 


and  1.0  within  the  desired  accuracy,  then  the  follow- 
ing replacement  is  made  to  restore  the  correct  normal- 
ization. 


This  normalization  check  and  replacement  strategy  con- 
tributes to  producing  a very  stable  and  accurate  numer- 
ical procedure. 


4.  Examples 

The  accuracy  and  power  of  the  developed  PIS  proce- 
dure will  be  illustrated  by  application  to  specific  case 
studies  taken  from  two  classes  of  dynamic  models.  Both 
models  are  described  by  Eq.  (3)  with  n=2  and  m=3. 
This  implies  a two-dimensional  state  space  vector 
X={XiJC2f^{Z^^.  Further,  TO (•)  and  Q(»)  are  such  that 
mi{Xi^2)=X2  and  q\j{»)=0  for y=l,2,3.  Assuming  suffi- 
cient restrictions  on  to(»)  and  0(*),  cf.  Refs.  [6,7],  X{t) 
becomes  a Markov  diffusion  process.  Invoking  Eq.  (13), 
it  can  be  shown  that,  up  to  correction  terms  of  order 
the  associated  TPD  assumes  the  form 


p{x,T  I x',0)=5(jCi— JC|'— JC2'  t)*P(x2,t  I x',0).  (25) 

p{x2,T  |x',0)  is  given  by  the  relation 


P(x2,t!x',0)= 


1 


V2tt/3(x')t 


exp 


(X2-X2 


'-TO2(x')t)^] 

2/3(x')t  y 


(26) 


Z+Zg{E}+h{Z)^rN{t).  (28) 

N{t)  denotes  a stationary,  zero-mean  Gaussian  white 
noise  satisfying  E[N(t}  iV(f+T)]=8(T),  where  8(»)  de- 
notes Dirac’s  delta  function,  Tis  a positive  constant  and 
g(£')  is  a function  of  the  total  energy  E=E(ZJZ)  given 
as  follows 

£=^Z'+V(Z)  (29) 

where 

Z 

Vizhj  h(s)ds.  (30) 

0 

For  this  example  m2(z,z)^-z  g[E(z,z)]—h(z),  and 

we  may  choose  g2i=q22=Wi=W2=0.  q23=r  and 

dW3(t)=N {t)dt . The  stationary,  joint  PDF,  denoted  by 
Pj(»),  is  then  determined  by  the  relation,  cf.  Refs. 
[20,21] 

E' 

Psiz,z)=C  exp{-  p j g(s)ds  },  (31) 

0 

where  £''=zV2+V(z),  and  C is  a normalization  constant 
to  ensure  a total  probability  equal  to  1.0. 

For  the  illustration  purposes  in  this  paper,  we  have 
chosen  the  following  special  case  of  Eq.  (28) 

Z(r)+2^(/){l+e[^ZV)+|z'(/)+  ^AZ\t)f^}+Z(t) 


where 

3 

^(x')=2  q2j(x')\  (27) 

j-i 

By  combining  Eqs.  (25)  and  (26),  and  applying  the 
solution  technique  described  in  the  previous  section,  the 
TPD  p(x,t  \x',t')  for  large  t—t'  can  be  calculated.  By 
this,  the  time  evolution  of  the  system  when  it  starts  from 
rest,  for  example,  can  be  studied.  The  stationary  PDF  is 
obtained  in  the  limit  as  r— For  application  of  the 
PIS  method  to  other  problems  involving  both  two-  and 
three-dimensional  state  space  vectors,  the  reader  may 
consult  Refs.  [11-13,18,19]. 

4.1  Example  1— The  Caughey  Oscillator 

There  is  a class  of  dynamic  models  for  which  there 
exist  an  analytical  solution  for  the  stationary  joint  PDF 
of  Z.  A member  of  this  class  may  be  called  a Caughey 
oscillator.  Ref.  [20],  The  generic  equation  of  motion  for 
this  oscillator  can  be  written  as 


-t-AZ'(r)=2V^(/)  (32) 

with  parameters  f,  e,  and  A. 

The  stationary  PDF  only  depends  on  the  parameters  e 
and  A , and  the  numerical  solution  for  the  following  set 
of  parameter  values  has  been  calculated  (e,A)=(0,0) 
(Gaussian  response),  (0,  0.2)  (Duffing  oscillator)  and 
(0.5,  0.1).  The  calculations  were  carried  out  with  the 
same  number  of  grid  points  on  both  axes  in  state  space, 
aviz.,  45.  Since  the  resulting  PDFs  are  actually  indepen- 
dent of  the  value  ^0. 1 was  chosen  for  the  Gaussian 
and  Duffing  cases,  while  ^0.5  was  adopted  for  the  last 
case.  The  time  increments  used  were  7M3.0025  s,  0.001 
s,  and  0.02  s,  respectively.  The  total  CPU  time  on  a DEC 
station  3100’  was  about  5 minutes  for  each  case.  In  Figs, 


'Certain  commercial  equipment,  instruments,  or  materials  are  identi- 
fied in  this  paper  to  specify  adequately  the  experimental  procedure. 
Such  identification  does  not  imply  recommendation  or  endorsement 
by  the  National  Institute  of  Standards  and  Technology,  nor  does  it 
imply  that  the  materials  or  'equipment  identified  are  necessarily  the 
best  available  for  the  purpose. 
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1 and  2 are  shown  the  marginal  PDFs  of  the  displace- 
ment response  for  the  three  case  studies  considered, 
together  with  the  corresponding  analytical  solutions.  In 
Fig.  3 are  given  the  corresponding  analytical  and  numer- 
ical results  for  the  mean  upcrossing  rate.  It  is  seen  that 
in  all  three  cases  the  agreement  between  the  numerical 
PIS  and  the  analytical  solution  is  very  good  over  the 
whole  range  of  probability  levels  given.  In  fact,  the 
accuracy  can  be  retained  down  to  much  lower  probabil- 
ity levels  (—10”"’)  at  a moderate  increase  in  computer 
time. 


are  independent  Gaussian  white  noises  satisfying 

E[Nj(t)Nj{t+T)]=rf  S(T),  >1,2,3,  (34) 

where  are  positive  constants.  For  this  example  it 
is  found  that  m2(z,z)=-2^z~y[z^+t^/(o^]-o)^z, 
q2\{z,z)=-2^zru  q22{z,z)=-oi^zr2  and  923(2 ,z)=T3. 

This  model  was  studied  by  Dimentberg  [22],  who 
showed  that  when 

oii  n (35) 


4.2  Example  2 — Parametric  and  External 
Excitation 


a closed-form  expression  for  the  stationary  joint  PDF 
can  be  obtained.  It  is  given  as 


In  this  example,  the  response  statistics  of  a nonlinear 
oscillator  subjected  to  both  external  and  parametric  ran- 
dom excitation  will  be  illustrated  by  applying  the 
methodology  of  the  paper  to  two  specific  case  studies. 

The  equation  of  motion  of  the  oscillator  is  the  follow- 
ing 


Z-h2^[  1 +N,{t )]Z+y [Z>  ^]Z+OA^[\ -HA^2(0]Z=>,(r ). 

<«>(i 

(33) 

Here  y,  and  tuu  are  positive  constants,  Nj{t),j=\, 2,3, 


Psl,Z,Z)  c > 


(36) 


where  C is  a normalization  constant  and 


a 


^ .4rQ  •>  ^ , ,2j^  “*■  2 ’ ^ .,.,2r'2 


(On 


(o^n 


(37) 


By  this,  we  have  the  opportunity  to  test  the  accuracy  of 
the  PIS  method  for  this  kind  of  dynamic  model.  The 
results  of  two  particular  cases  will  be  presented. 


Displacement 

Fig.  1.  FYobability  density  function  of  displacement  response  for  the  Caughey  oscillator  in  example  1.  Ana- 
lytical solutions:  , e=0,  A=0;  - • - e=0,  A=0.2; , £=0.5,  A=0.l.  Numerical  path  integral 

solution:  -1-,  £=0,  A=0;  X,  e=0,  A=0.2;  O,  £=0.5,  A=0. 1. 
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Fig.  2.  Logarithmic  plot  of  the  probability  density  function  of  displacement  response  for  the  Caughey  oscillator 
in  example  1.  Key  as  in  Fig.  1. 


Displacement 


Fig.  3.  Mean  upcrossing  rate  of  displacement  response  for  the  Caughey  oscillator  in  example  1 . Key  as  in  Fig. 
1. 
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Case  1:  Here  the  following  parameter  values  were 
used.  ^0.1,  y=0.1,  <u»,=  1.0,  /7=2.5,  71=0.1,  71=0.3. 
For  the  numerical  calculations  a grid  size  of  49X49 
points  and  a time  increment  t=0.01  s was  used.  The  total 
CPU  time  on  a DEC  3100  work  station  was  3 min  for  the 
PIS  calculation.  The  results  for  the  analytical  and  nu- 
merical solutions  are  given  in  Figs.  4-6.  In  Figs.  4 and 
5 are  shown  the  marginal  PDF  of  the  displacement  re- 
sponse and  in  Fig.  6 is  shown  the  corresponding  mean 
upcrossing  rate. 

Case  2:  In  this  case  the  following  set  of  parameters 
were  used.  ^0.1,  y=0.4,  <uc)=1.0,  71=5.0,  71=0.2, 
71=0.3.  A grid  size  of  51X51  points  together  with  a 
time  increment  t=0.01  s were  chosen.  The  CPU  time 
was  about  the  same  as  in  the  previous  case.  The  same 
results  as  for  Case  I are  presented  in  Figs.  4-6. 


5.  Conclusions 

A numerical  method  for  estimating  the  extreme  re- 
sponse of  nonlinear  oscillators  excited  by  white  noise, 
or  filtered  white  noise,  has  been  described.  The  example 
calculations  presented  show  that  the  method  gives  very 
accurate  estimates  of  the  required  joint  PDF.  In  fact,  for 
every  example  having  analytical  solution  on  which  the 
method  has  been  tested,  complete  agreement  has  been 
found  with  proper  choice  of  grid  size  and  time  incre- 
ment in  the  numerical  solution  procedure.  In  the  present 
paper,  of  course,  only  a few  cases  can  be  given.  Experi- 
ence with  the  method  indicates  that  two-dimensional 
problems  can  be  solved  routinely  with  high  accuracy 
requiring  a few  minutes  CPU  time  on  a work  station 
(DEC  station  31(X)).  The  solution  of  three-dimensional 
problems  requires  more  care  in  the  sense  that  computer 
capacity  starts  to  become  an  issue  of  importance.  In 
such  cases  the  CPU  time  easily  runs  into  hours. 


Fig.  4.  Probability  density  function  of  displacement  response  for  the  oscillator  in  example  2,  case  1 and  2. 
Analytical  solutions:  , case  case  2.  Numerical  path  integral  solution:  +,  case  1;  O,  case  2. 
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Fig.  5.  Logarithmic  plot  of  the  probability  density  function  of  displacement  response  for  the  oscillator  in 
example  2,  case  1 and  2.  Key  as  in  Fig.  4. 


Fig.  6.  Mean  upcrossing  rate  of  displacement  response  for  the  oscillator  in  example  2,  case  1 and  2.  Key  as  in 
Fig.  4. 
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